%Simulate the future state of each hospital by starting from the actual
%state each year, including decreasing sunk cost.


tic

choice=1:13;
ns=numel(choice);
beta=0.95;
tbar=3;
%The alpha and gamma are the coefficients from the policy function for
%endogenous hospitals (stand-alone hospitals) that haven't adopted EMR.
%alpha=[0.0163	0.665	0	0	0.476]; %the coefficients from the policy function (vary across choices and hospitals): vmkt vmsh vlag vsys vmkt*big

alpha=[0.0888	0.830	0.983	1.038 0 0 0 0  0.0558];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q vdom*yrsince

mktcat_alpha=[0.972	1.888	4.209];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_alpha=[zeros(numel(mktcat_alpha),1),repmat(mktcat_alpha',1,ns-1)];

gamma2=[0.000572	-0.757	-15.50	0.266	0.128	0.826	-9.734]; %the following gamma's are coefficients for the choice-fixed-hospital-varied variables (from estimating the policy function based on the entire sample)
gamma3=[0.00243	0.354	0.857	2.522	0.516	-0.111	-12.12];  %the order is as follows: preadp(indicate whether adopting or not) bdtot nprofit teaching HHI totpop65 yrsince constant 
gamma4=[-0.0110	0.0433	-15.24	2.049	0.228	0.163	-6.617];
gamma5=[-0.0113	0.0279	-15.20	1.014	-0.116	0.0254	-2.953];
gamma6=[0.00874	1.331	-16.41	-1.867	0.0276	-0.329	-7.629];
gamma7=[0.00304	1.046	1.825	-0.0634	-0.0397	-0.443	-6.071];
gamma8=[0.00933	1.621	-0.901	-1.019	-0.400	-0.589	-3.885];
gamma9=[-0.00635	-0.423	-15.22	0.342	0.0396	0.247	-5.481];
gamma10=[0.00513	0.280	-1.811	2.640	0.269	0.0362	-9.106];
gamma11=[0.00506	1.172	-1.472	1.614	0.371	-0.318	-10.22];
gamma12=[0.000176	0.454	-1.850	2.256	0.311	0.323	-8.721];
gamma13=[0.00120	0.493	0.118	2.467	0.293	0.113	-8.938];


Gamma=[zeros(size(gamma2))' gamma2' gamma3' gamma4' gamma5' gamma6' gamma7' gamma8' gamma9' gamma10' gamma11' gamma12' gamma13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)


%The coefficients for stand-alone hospitals that have adopted EMR.
tau=[-0.0719	0.630	1.262	0.565	6.221	6.235	6.009	5.714	-0.0429]; %the coefficients from the policy function (vary across choices and hospitals): vmkt vmsh vlag vsys vlag*big vdom*yrsince
mktcat_tau=[0	0	0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

zeta2=[0 0 0 0 0 0 0]; %the following gamma's are coefficients for the choice-fixed-hospital-varied variables (from estimating the policy function based on the entire sample)
zeta3=[0.00110	-1.312	0.264	-0.0646	0.0467	-0.0632	0.564 ];  %the order is as follows: preadp(indicate whether adopting or not) bdtot profit teaching HHI totpop65 yrsince constant 
zeta4=[-0.0121	0.986	-10.64	1.587	0.0892	-0.427	1.755];
zeta5=[-0.0199	-1.061	-10.80	0.687	-0.0455	0.164	3.141];
zeta6=[0.000817	-0.896	1.162	2.637	0.334	0.152	-4.395];
zeta7=[ 0.00150	-2.820	0.952	1.106	0.189	-0.555	-0.563];
zeta8=[0.00189	-0.132	1.279	4.625	0.873	0.415	-13.79];
zeta9=[-0.0108	-0.977	-10.50	-1.054	-0.297	-0.0195	5.397];
zeta10=[ 0.000797	-2.284	0.268	-0.620	-0.306	-0.432	5.185];
zeta11=[0.00124	-0.616	0.126	0.452	0.0902	0.0492	-1.019];
zeta12=[-0.00117	-2.174	-1.460	1.107	0.163	0.0640	0.229];
zeta13=[-0.000625	0.0622	0.609	0.487	0.0981	-0.188	-1.311];

Zeta=[zeros(size(zeta2))' zeta2' zeta3' zeta4' zeta5' zeta6' zeta7' zeta8' zeta9' zeta10' zeta11' zeta12' zeta13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)



%The sigma and eta are the coefficients from the policy function based on
%the sample of exogenous hospitals (affiliated hospitals) that haven't adopted EMR.
sigma=[0 0 0 2.843]; %the coefficient for vsysh 

mktcat_sigma=[0.816	1.724	4.466];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_sigma=[zeros(numel(mktcat_sigma),1),repmat(mktcat_sigma',1,ns-1)];

eta2=[0.00116	-0.725	-0.522	2.598	0.505	0.623	-12.97]; %the order is as follows: preadp(indicate whether adopting or not) bdtot nprofit perc_mcr hhi totpop65 yrsince constant
eta3=[0.00191	0.283	-1.152	0.506	0.123	-0.356	-4.910];
eta4=[-0.00903	-0.607	2.212	0.451	0.227	0.0796	-7.879];
eta5=[-0.0128	-0.265	0.541	1.338	-0.0408	-0.0633	-4.130];
eta6=[0.00143	1.816	-2.623	4.811	0.983	-0.557	-17.57];
eta7=[0.00144	3.075	-1.358	1.166	0.196	-0.808	-7.753];
eta8=[0.00237	2.194	-0.720	5.653	0.818	-0.560	-16.97];
eta9=[-0.00216	-1.185	-2.916	1.476	-0.0537	-0.531	-1.384];
eta10=[0.00190	0.514	-1.433	1.819	0.328	-0.329	-7.590];
eta11=[0.000290	0.0908	-3.575	0.355	0.129	-0.293	-3.749];
eta12=[-0.0000573	1.084	-1.179	1.318	0.133	-0.177	-5.849];
eta13=[0.00146	0.824	-2.148	0.558	-0.0441	-0.301	-3.594];

Eta=[zeros(size(eta2))' eta2' eta3' eta4' eta5' eta6' eta7' eta8' eta9' eta10' eta11' eta12' eta13'];



%The coefficients for the affiliated hospitals that have adopted EMR.
rho=[0 0 4.404   2.244]; %the coefficient: vmkt vmsh vlag vsys

mktcat_rho=[0    0   0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

lambda2=[0 0 0 0 0 0 0]; %the order is as follows: preadp(indicate whether adopting or not) bdtot perc_mcr perc_mcd hhi totpop65 yrsince constant

lambda3=[0.00175	4.419	3.255	-1.351	-0.122	0.296	-1.051];
lambda4=[-0.00511	4.138	3.765	-2.602	-0.595	0.508	4.475];
lambda5=[-0.00634	4.813	1.974	-5.534	-1.125	0.346	10.43];
lambda6=[0.00228	3.999	1.493	-1.266	0.108	0.349	-4.069];
lambda7=[0.00274	3.788	-0.828	-3.348	-0.336	0.263	2.884];
lambda8=[0.00116	2.588	0.123	-1.254	-0.0860	0.523	-1.334];
lambda9=[-0.00463	4.149	4.885	-1.418	-0.315	0.212	1.266];
lambda10=[0.000645	4.501	5.371	-2.815	-0.312	0.334	0.764];
lambda11=[-0.0000678	1.360	2.449	-0.495	-0.416	0.589	2.338];
lambda12=[-0.000544	4.215	3.890	-1.108	-0.261	0.620	0.216];
lambda13=[0.000438	7.017	5.461	-0.993	-0.106	0.950	-4.787];
Lambda=[zeros(size(lambda2))' lambda2' lambda3' lambda4' lambda5' lambda6' lambda7' lambda8' lambda9' lambda10' lambda11' lambda12' lambda13'];



bed25=25;
bed50=80;
bed75=193;
%ns=1;       %will take 516 secs if ns=500
T=100;
%s=1;
textFileName3=sprintf('fs%d_decsunk.csv',s);  %the output file
textFileName4=sprintf('ddom%d_decsunk.csv',s);   %the output file
textFileName5=sprintf('eerr%d_decsunk.csv',s);   %the output file
textFileName6=sprintf('mks%d_decsunk.csv',s);   %the output file

rng(s);
FSTA=[];
DDOM=[];
EERR=[]; %### UPDATED on 03/16/2017: the output including the probability as the expected error
MKTS=[]; %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation

textFileName_actmkt=sprintf('mktvar_decsunk.csv');  %the input file recording actual market variables
actmkt=dlmread(textFileName_actmkt);  %the matrix: ahaid year vid vmshbd* vlagbd* vdom_yrsince newvsysh vlag hsa sysid


for k=6:9
%k=7;
textFileName1a=sprintf('prdprob_a%d_decsunk.csv',k); %the input file
textFileName1na=sprintf('prdprob_na%d_decsunk.csv',k); %the input file
textFileName2a=sprintf('syschs_a%d_decsunk.csv',k);  %the input file
textFileName2na=sprintf('syschs_na%d_decsunk.csv',k);  %the input file
J1=dlmread(textFileName1na);  % ### UPDATED on 04/21/2020---the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName1a); 
%reshape the raw data from STATA, SA-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
H1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, SA-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
H2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
P=[H1;H2];
%pick=find(P(:,2)==17 | P(:,2)==10 |P(:,2)==44 |P(:,2)==46 |P(:,2)==62 |P(:,2)==90);
%P=P(pick,:);
J1=dlmread(textFileName2na); %the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName2a); 
%reshape the raw data from STATA, AF-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
HH1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, AF-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
HH2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
PP=[HH1;HH2];
allsys=PP(:,end);
%Specify each page of FSTA and DDOM: fs and G
fs=zeros(size(P,1)*size(choice,2),T+5); 
G=zeros(size(P,1)*size(choice,2),T+4);
MS=zeros(size(P,1)*size(choice,2),T+4); %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation
eer=zeros(size(P,1)*ns,T+3);  %### UPDATED on 03/16/2017: include the simulated expected errors: ahaid, year, option given, prob@t=1, prob@t=2, ..., prob@t=100

for i=1:size(P,1)  
    tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
    exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
    exomkt=sortrows(exomkt,1);
    altemp=[exomkt;tempHH];  %group all the rival hospitals (standalone and affiliated) together (to do the random draw)
    tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
    alltemp=[tempmkt;tempHH]; %group all the hospitals (including the one considered) (standalone and affiliated) together (to do the random draw)
    %The following is to compute fs0 and g0. Then add [year, bdtot, fs0]
    %into fs and add [year,g0] into G.
    fs0=P(P(:,1)==P(i,1),4);
    mkts0=0;
    prechs=altemp(:,4);
    prechs(prechs==1)=[];
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
       mktsh=occur/sum(occur);
       mkts0=mktsh(fs0);
    end    
    prechs(prechs==2)=[];
    prechs(prechs==13)=[];
    if isempty(prechs)==1
       g0=0;
    end
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
     %  if occur(2)>0
     %     occur(2)=1;
     %  end
       occur=reshape(occur,1,[]);
       mktdom=choice(occur==max(occur));
       g0=max(fs0-mktdom==0);
    end
  %%%%%%%%%%%%%%%%%%%%%%%%%
  %consider the case where the hospital has already installed EMR  
   if P(i,4)>1   
        want=zeros(size(choice,2)-1,T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2)-1,T+1);  
        mkts=zeros(size(choice,2)-1,T+1);
        tempeer=zeros(ns-1,T+1);  %### UPDATED on 03/16/2017: a section of eer: the probability at each period
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market. ### 04/21/2020: should be
        %for stand-alone hospitals whose competitors are ALL stand-alone
        %hospitals, too.
        if isempty(sys)==1
          for l=2:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;
              if isempty(altemp)==1
                  newchs=l;
              end
              if isempty(altemp)==0
                 %UPDATED on 04/26/2020: calculate the probability of each choice at the base year (i.e., the actual year)    
                 actv1=zeros(size(altemp,1),ns);
                 actv2=zeros(size(altemp,1),ns);
                 for ii=1:size(altemp,1)
                     regressor=actmkt(actmkt(:,1)==altemp(ii,1)&actmkt(:,2)==2000+k,4:12);
                     actv1(ii,:)=(regressor*alpha')'+altemp(ii,[3,18,20,21,22])*Gamma([1,2,3,4,5],:)...
                         +altemp(ii,[25,26,27])*mktcat_alpha...
                         +tbar*Gamma(6,:)+Gamma(end,:);                     
                     actv2(ii,:)=(regressor*tau')'+altemp(ii,[3,19,20,21,22])*Zeta([1,2,3,4,5],:)...
                         +tbar*Zeta(6,:)+Zeta(end,:);                     
                 end                
                 ind1_actv=(altemp(:,4)==1);
                 ind2_actv=(altemp(:,4)>1);
                 %ind3=(univ(:,end)>0&univ(:,4)==1);
                 %ind4=(univ(:,end)>0&univ(:,4)>1);
                 actv=actv1.*ind1_actv(:,ones(size(actv1,2),1))+actv2.*ind2_actv(:,ones(size(actv2,2),1));
                 %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                 exp_actv=exp(actv);
                 exp_actv(altemp(:,4)>1,1)=0;
                 sumexp_actv=sum(exp_actv,2);
                 dom_expactv=sumexp_actv(:,ones(size(exp_actv,2),1));
                 w=exp_actv./dom_expactv;                  
                 %w=altemp(:,5:17);
                 cc=cumsum(w,2);
                 cc(:,end)=1;
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1  %for monopoly market
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l;
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                                        
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0;
                       %mkts(l-1,t+1)=0;
                     end
                     if isempty(prechs)==0
                        occur=histc(prechs,choice)';
                        occur=reshape(occur,1,[]);           
                        vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                        %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                        %vmshf=vmshff;
                        prechs(prechs==2)=[];
                        prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0; %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4);
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);  %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".                   
                     Tchoice=choice';                     
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);
                    
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';  
                     % ### UPDATED on 04/26/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:)...
                       +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:)...
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);      
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                                        
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);  %###!!!If the first one changed back to using max, remember to use "altemp"!
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];                   
                     %### UPDATED on 1/4/2017: re-define "selfnewchs" by
                     %figuring out what brings in largest payoff given the
                     %rivals' new choice, from "newchs(2:end)".            
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),2000+k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),(2000+k)*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market. ### 04/21/2020: consider hospitals whose competitors
        %include affiliated hospitals
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:); %find out the outside ones
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=2:13
              %w=altemp(:,5:17);
              %chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);              
              %UPDATED on 04/26/2020: calculate the probability of each choice at the base year (i.e., the actual year)    
              actv1=zeros(size(exouniv,1),ns);
              actv2=zeros(size(exouniv,1),ns);
              actv3=zeros(size(exouniv,1),ns);
              actv4=zeros(size(exouniv,1),ns);              
              for ii=1:size(exouniv,1)
                  regressor=actmkt(actmkt(:,1)==exouniv(ii,1)&actmkt(:,2)==2000+k,4:12);
                  regressor_sys=actmkt(actmkt(:,1)==exouniv(ii,1)&actmkt(:,2)==2000+k,13:14);
                  actv1(ii,:)=(regressor*alpha')'+exouniv(ii,[3,18,20,21,22])*Gamma([1,2,3,4,5],:)...
                         +exouniv(ii,[25,26,27])*mktcat_alpha...
                         +tbar*Gamma(6,:)+Gamma(end,:);                     
                  actv2(ii,:)=(regressor*tau')'+exouniv(ii,[3,19,20,21,22])*Zeta([1,2,3,4,5],:)...
                         +tbar*Zeta(6,:)+Zeta(end,:);                      
                  actv3(ii,:)=(regressor_sys*sigma(end-1:end)')'+exouniv(ii,[3,18,23,21,22])*Eta([1,2,3,4,5],:)...
                         +exouniv(ii,[25,26,27])*mktcat_sigma...
                         +tbar*Eta(6,:)+Eta(end,:);                                                
                  actv4(ii,:)=(regressor_sys*rho(end-1:end)')'+exouniv(ii,[3,23,24,21,22])*Lambda([1,2,3,4,5],:)...
                         +tbar*Lambda(6,:)+Lambda(end,:);                                    
              end             
              ind1_actv=(exouniv(:,end)==0&exouniv(:,4)==1);
              ind2_actv=(exouniv(:,end)==0&exouniv(:,4)>1);
              ind3_actv=(exouniv(:,end)>0&exouniv(:,4)==1);
              ind4_actv=(exouniv(:,end)>0&exouniv(:,4)>1);
              actv=actv1.*ind1_actv(:,ones(size(actv1,2),1))+actv2.*ind2_actv(:,ones(size(actv2,2),1))+...
                actv3.*ind3_actv(:,ones(size(actv3,2),1))+actv4.*ind4_actv(:,ones(size(actv4,2),1));
              exp_actv=exp(actv);
              exp_actv(exouniv(:,4)>1,1)=0;
              sumexp_actv=sum(exp_actv,2);
              dom_expactv=sumexp_actv(:,ones(size(exp_actv,2),1));
              w=exp_actv./dom_expactv;         
              %w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1;
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l; 
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     
                    alltemp(:,4)=newchs(1:size(alltemp,1),:);
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                     
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0; %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0;  %### UPDATED on 07/16/2018                     
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                  
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                    vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                    vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                    

                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);                           
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);                    
                     %bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);                     
                     %bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)>1);                                      
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                 
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                     % +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                     % +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                        
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                    

                     % ### UPDATED on 04/26/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:)...
                       +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:)...                  
                       +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:);                                       
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:)...                                   
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);                          
                     v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:)...
                      +yrsince*Eta(6,:);                                                           
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:)...
                      +yrsince*Lambda(6,:);
                                        
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                     
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);  %!!! change it to "univ" or "exouniv" depends on random or max.
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),2000+k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),(2000+k)*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
   end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %consider the case where the hospital has never installed EMR
   if P(i,4)==1  
        %tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
        %exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
        %exomkt=sortrows(exomkt,1);
        %altemp=[exomkt;tempHH];  %group all the hospitals (standalone and affiliated) together to do the random draw)
        %tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pc1 pc2...pc13
        %alltemp=[tempmkt;tempHH];
        want=zeros(size(choice,2),T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2),T+1); 
        mkts=zeros(size(choice,2),T+1);
        tempeer=zeros(ns,T+1); 
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market.
        if isempty(sys)==1
          for l=1:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;             
              if isempty(altemp)==1  %consider the case of monopoly
                  newchs=l;
              end
              if isempty(altemp)==0
                 %UPDATED on 04/26/2020: calculate the probability of each choice at the base year (i.e., the actual year)    
                 actv1=zeros(size(altemp,1),ns);
                 actv2=zeros(size(altemp,1),ns);
                 for ii=1:size(altemp,1)
                     regressor=actmkt(actmkt(:,1)==altemp(ii,1)&actmkt(:,2)==2000+k,4:12);
                     actv1(ii,:)=(regressor*alpha')'+altemp(ii,[3,18,20,21,22])*Gamma([1,2,3,4,5],:)...
                         +altemp(ii,[25,26,27])*mktcat_alpha...
                         +tbar*Gamma(6,:)+Gamma(end,:);                     
                     actv2(ii,:)=(regressor*tau')'+altemp(ii,[3,19,20,21,22])*Zeta([1,2,3,4,5],:)...
                         +tbar*Zeta(6,:)+Zeta(end,:);                     
                 end                
                 ind1_actv=(altemp(:,4)==1);
                 ind2_actv=(altemp(:,4)>1);
                 %ind3=(univ(:,end)>0&univ(:,4)==1);
                 %ind4=(univ(:,end)>0&univ(:,4)>1);
                 actv=actv1.*ind1_actv(:,ones(size(actv1,2),1))+actv2.*ind2_actv(:,ones(size(actv2,2),1));
                 %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                 exp_actv=exp(actv);
                 exp_actv(altemp(:,4)>1,1)=0;
                 sumexp_actv=sum(exp_actv,2);
                 dom_expactv=sumexp_actv(:,ones(size(exp_actv,2),1));
                 w=exp_actv./dom_expactv;                  
                 %w=altemp(:,5:17);
                 cc=cumsum(w,2);
                 cc(:,end)=1;                 
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];
                                          
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;   %### UPDATED on 07/16/2018
                     end
                     
                     
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);     %### UPDATED on 07/16/2018                 
                     Tchoice=choice';      
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                   
                   curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);


                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);    
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);                  
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);                  
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                     
                     % ### UPDATED on 04/26/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:)...
                       +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:)...
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);      
                      %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;              
                     %selfnewchs=find(w(1,:)==max(w(1,:)));
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l,t+1)=newchs(1);  
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]; 
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];  %### UPDATED on 07/16/2018
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),(2000+k)*ones(size(tempeer,1),1),tempeer];
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market.
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that have belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:);
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=1:13
              %w=altemp(:,5:17);
              %chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              
              %UPDATED on 04/26/2020: calculate the probability of each choice at the base year (i.e., the actual year)    
              actv1=zeros(size(exouniv,1),ns);
              actv2=zeros(size(exouniv,1),ns);
              actv3=zeros(size(exouniv,1),ns);
              actv4=zeros(size(exouniv,1),ns);              
              for ii=1:size(exouniv,1)
                  regressor=actmkt(actmkt(:,1)==exouniv(ii,1)&actmkt(:,2)==2000+k,4:12);
                  regressor_sys=actmkt(actmkt(:,1)==exouniv(ii,1)&actmkt(:,2)==2000+k,13:14);
                  actv1(ii,:)=(regressor*alpha')'+exouniv(ii,[3,18,20,21,22])*Gamma([1,2,3,4,5],:)...
                         +exouniv(ii,[25,26,27])*mktcat_alpha...
                         +tbar*Gamma(6,:)+Gamma(end,:);                     
                  actv2(ii,:)=(regressor*tau')'+exouniv(ii,[3,19,20,21,22])*Zeta([1,2,3,4,5],:)...
                         +tbar*Zeta(6,:)+Zeta(end,:);                      
                  actv3(ii,:)=(regressor_sys*sigma(end-1:end)')'+exouniv(ii,[3,18,23,21,22])*Eta([1,2,3,4,5],:)...
                         +exouniv(ii,[25,26,27])*mktcat_sigma...
                         +tbar*Eta(6,:)+Eta(end,:);                                                
                  actv4(ii,:)=(regressor_sys*rho(end-1:end)')'+exouniv(ii,[3,23,24,21,22])*Lambda([1,2,3,4,5],:)...
                         +tbar*Lambda(6,:)+Lambda(end,:);                                    
              end             
              ind1_actv=(exouniv(:,end)==0&exouniv(:,4)==1);
              ind2_actv=(exouniv(:,end)==0&exouniv(:,4)>1);
              ind3_actv=(exouniv(:,end)>0&exouniv(:,4)==1);
              ind4_actv=(exouniv(:,end)>0&exouniv(:,4)>1);
              actv=actv1.*ind1_actv(:,ones(size(actv1,2),1))+actv2.*ind2_actv(:,ones(size(actv2,2),1))+...
                actv3.*ind3_actv(:,ones(size(actv3,2),1))+actv4.*ind4_actv(:,ones(size(actv4,2),1));
              exp_actv=exp(actv);
              exp_actv(exouniv(:,4)>1,1)=0;
              sumexp_actv=sum(exp_actv,2);
              dom_expactv=sumexp_actv(:,ones(size(exp_actv,2),1));
              w=exp_actv./dom_expactv;         
              %w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1; 
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     alltemp(:,4)=newchs(1:size(alltemp,1),:);
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];                               
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;   %### UPDATED on 07/16/2018                   
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                 
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                     vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                     vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                   
                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);  
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);    
                                     
                    % bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);     
                    % bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)==1);                                        
                    % v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                    %  +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                   
                    % v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                    %  +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                    % v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                    %  +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                    % v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                    %  +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                    
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';
 
                     % ### UPDATED on 04/26/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:)...
                       +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:)...                  
                       +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:);                                       
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:)...                                   
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);                          
                     v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:)...
                      +yrsince*Eta(6,:);                                                           
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:)...
                      +yrsince*Lambda(6,:);
                                        
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;
                     %{
                     selfnewchs=find(w(1,:)==max(w(1,:)));
                     w(1,:)=[];
                     cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l,t+1)=newchs(1); 
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g];   
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),(2000+k)*ones(size(tempeer,1),1),tempeer];

        end
    end
end
G(:,[3,4])=G(:,[4,3]);  %exchange g0 with action_given
MS(:,[3,4])=MS(:,[4,3]); % ### UPDATED on 08/22/2018: exchange mkts0 with action_given
FSTA=[FSTA;fs];
DDOM=[DDOM;G];
MKTS=[MKTS;MS];
EERR=[EERR;eer];
end
dlmwrite(textFileName3,FSTA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName4,DDOM,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName5,EERR,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName6,MKTS,'delimiter', ',', 'precision', 8)


toc


